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Abstract. A three-dimensional model of the far turbulent wake behind a self-propelled body 
in a passively stratified medium is considered. The model is reduced to a system of ordinary 
differential equations by a similarity reduction and the i?-determining equations method. 
The system of ordinary differential equations satisfying natural boundary conditions is solved 
numerically. The solutions obtained here are in close agreement with experimental data. 
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1 Introduction 

Most flows occurring in nature and engineering practice are turbulent (see, e.g., [10, 20, 22]). 
Semiempirical models of turbulence are widely used in the modeling of turbulent flows 
[16, 21, 27]. However, there are only a few analytical results on such models (see, e.g., [2, 11]). 

The far turbulent wake behind an axisymmetric body in a stratified medium is an example of 
a free shear flow. Sufficiently complete experimental data on the dynamics of turbulent wakes 
generated by moving bodies in stratified fluids were obtained by Lin and Pao [17]. 

The far turbulent wake behind an axisymmetric towed body in a linearly stratified medium 
was numerically simulated in [9]. Chernykh et al. [6] carried out the numerical simulation of the 
dynamics of turbulent wakes in a stable stratified medium based on hierarchy of second order 
closure models. Calculation results obtained in [6, 9] are in close agreement with experimental 
data [17]. 

Similarity solutions for several turbulence models were constructed in [7, 13, 14, 15]. In 
the current paper we consider three-dimensional semiempirical model of the far turbulent wake 
behind an axisymmetrical self-propelled body in a passively stratified medium (see [6, 4, 25] 
and the references therein). 

This paper is organized as follows. In Section 3 we determine the most general continuous 
classical symmetry group of the model and obtain the similarity reduction of the model. In 
Section 4 we use the 5-determining equations (BDEs) method [1, 12] to transform the reduced 
system into a system of ordinary differential equations (ODEs). 

In the last section, we consider a boundary value problem for the system of ODEs. We use 
the modified shooting method and the asymptotic expansion of the solution in the vicinity of 
the singular point to solve this problem. Finally, computational results are given. 



*This paper is a contribution to the Special Issue "Geometrical Methods in Mathematical Physics". The full 
collection is available at http://www.emis.de/journals/SIGMA/GMMP2012.html 
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2 Model 

The following three-dimensional semiempirical model of turbulence was constructed in [4, 6, 25] 
to calculate characteristics of the far turbulent wake behind an axisymmetric self-propelled body 
in a passively stratified medium: 

^ 9e d ^ de d ^ de 

ox oy e oy oz e oz 

Bf B Bf B p*^ Bf 

^0?^ = tCe-^ + tCe-^ - Ce2-, (2) 

Bx By e By Bz e Bz e 

Uq^^ = —C —"^^ + —C —"^^ - —C — (3) 
Bx By ^ e By Bz ^ e Bz Bz ^ e ^ 
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Bx By e By Bz e Bz e By 

where e is the turbulent kinetic energy, e is the kinetic energy dissipation rate, {pi) is the averaged 
density defect, and (/?'^) is the density fluctuation variance. All the unknown functions depend 
on X, y, and z. The quantities Ce = 0.136, Q = Ce/S, 6 = 1.3, Ce2 = 1-92, Cp = 0.208, 
Cip = 0.087, Ct = 1-25 are generally accepted empirical constant [8, 21]. Uq is the free stream 
velocity. The marching variable x in the equations (1)~(4) acts as time. 

This model is based on the three-dimensional parabolized system of averaged Navier-Stokes 
equations in the Oberbeck-Boussinesq approximation (see [5, 26]) 

BUd BUd , ^-.rdUd B , , B , , 

BV ^,,dV BV 1 B{pi) a 2. d 
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where Ud = Uq — U is the defect of the averaged longitudinal velocity component; U, V and W 
are the mean flow velocity component along x-, y- and z-axes, respectively; {pi) is the deviation 
from the hydrostatic pressure due to stratification Ps{z); g is the gravity acceleration; (pi) is 
the averaged density defect: pi = p — ps', Ps = Ps{z) is the undisturbed fluid density assumed to 
be linear: Ps{z) = po{l — az), a > is a constant; the prime indicates the pulsating components; 
( ) indicates averaging. 

In [9, 26] the Reynolds stress tensor components {u[u'j), the turbulent flows {u[p'), and the 
density fluctuation variance (p'^) are defined by the algebraic relations [21]. Since the flow in 
the far turbulent wake is considered, these relations are simplified as follows 
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Differential transport equations [21] are used in [4, 6, 25] to determine the kinetic turbulent 
energy e, the kinetic energy dissipation rate e, and the shear Reynolds stress v'w': 



Uo^ + + = + ^Kj/ + P + G-e, (20) 
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\ dy dz J po e 

The turbulent viscosity coefficients in these equations are K^y = Ky, K^z = Kz, Kgry = K^y/S, 
Ksz = Kez/S. The model (l)-(4) is an analogue of the equations (5)-(22) (for the diffusion 
approximation in a homogeneous fluid ^ = 0, 1^ = 0, and g = 0). In what follows, we assume 
that the free stream velocity Uq equals unity. 



3 Similarity solution 

The infinitesimal symmetry group [18, 19] of the model (l)-(4) is spanned by the eight vector 
fields 

d d d d d d d 
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Available experimental data [10, 17] and numerical calculations [9, 20, 27] show that the flow in 
the far turbulent wake can be considered to be close to a self-similar flow. We therefore consider 
the linear combination of scaling vector fields Xq and Xj 

The similarity variables associated with the infinitesimal generator Z are 

where G, H, and R are arbitrary functions. According to physical considerations (the 
influence of gravity is neglected in this problem) , the functions E and G must be of the form 

E{tr]) = E(Ve + V^), G{C,v) = G(Ve + V^)- (23) 



We obtain the reduced system by introducing similarity variables. Using (23) and changing 
to polar coordinates ^ = r cos (j), r] = r sin (p the reduced system becomes 



Ce^ {eE" + 2E'^ - ^E'G' + ^E'^ + arE' + 2(1 -a)E-G = 0, 
Ce^ (eg" - + 2E'G' + -G'] + arG' + (3 - 2a)G - C,2^ = 0, 

Cp^ (h„ + ^H.^i + (g.^ (2E' -^G' + -]+ar] K,. - aH 



(24) 
(25) 



+ Cp- ( -G'-2^' ) sin0 = O, (26) 



G \ ''"^ ^2 <pv I ^ x^G \ G r I 

E fE 
'G [g 

Cip^ [Rrr + ^i?<^<^) + (Ci,^ + 2E' - ^G') + ar) R^. - [ct^ + 2a) R 

+ 2Cp^ {h^, + ^Hl^ + 2Cp^ - AGp^ (^Hr sin </> + H^^^ = 0, (27) 

where E = E{r), G = G{r), H = H{r,(j)), and R = R{r,(j)). Here, and throughout, subscripts 
denote derivatives, so Hr = dH/dr, etc. 

Lie's classical method do not provide solution of the reduced system agreed with experimental 
data. We therefore use the BDEs method. 

4 BDEs method 

The concept of BDEs of a system of partial differential equations (PDEs) was introduced 
in [1, 12]. Consider a scalar PDE 

^7(x,^i,n^^t^...) = 0, (28) 

where x = {xi,X2, ■ ■ ■ ,Xn) denotes n independent variables, u denotes the dependent variable, 
and 

Uk 
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denotes the set of coordinates corresponding to all k-th order partial derivatives of u with respect 
to X. In the BDEs method, an extension of the classical symmetry determining relations is made 
by incorporating an additional factor i?(x, u, n^, n^, . . . ). For a scalar PDE (28), BDE is 

+ E ^°(^)|^ + Bh\n=o = 0, (29) 

where /i is a function of x, u, v},v?', . . . ; = D"^^ ■ ■ ■ D^^, Dx^ is the total Xi derivative; a is 
a multi-index. Equality (29) must hold for all solutions of (28). 

Now we use the BDEs method to reduce (26) and (27) to some ODEs. Consider more general 
equation than (26) 

H^^ + r^Hrr + A{r)Hr + B{r)H + C{r) sin (/> = 0, (30) 

where A{r), B{r), and C(r) are arbitrary functions. BDE corresponding to (30) is 

Dlh + r'^Djh + 6i(r, (/.)L>^/i + b2{r, (l))h = 0. (31) 

Here, and throughout Z)^, Dr are the operators of total differentiation with respect to and r. 
The function h may depend on r, (p, H and derivatives of H. The functions bi{r, (p) and b2{r, 0) 
are to be determined together with the function h. Note that if we let in (31) 

bi{r,(l)) = A{r), b2{r,(t>) = B{r), 

we obtain the classical determining equation [18, 19]. 

For simplicity, we assume that a solution of (31) is independent of r and partial derivatives 
of H with respect to r 

h = H^^ + hi (0, H, H^) . (32) 

Substituting (32) into (31) gives a polynomial equation for derivatives of the fourth order. This 
polynomial must identically vanish. We can express the derivatives Hrr(f,cf>, H^^^^, Hr(j,(j,, H^^^^ 
and using (30). The coefficient of Hrrr implies bi{r,(j)) = A{r). 

As a result, the left side of (31) is a polynomial in Hrr and Hnj). Collecting similar terms we 
obtain 

2 {A{r)Hr + B{r)H + C(r) sin 0) h^^^^ - 2H^hi^^, - 2h^„^ + B{r) - b^ir, 0) = 0, 
hiH^H^ = 0, hi^^^ = 0. 
Hence 

hi (0, H, H^) = h2{(l))H^ + hs{^, H), b2{r, 0) = B{r) - 2h'^{^). 

Substituting the functions hi, b2, and hi into the left side of (31) we obtain the polynomial 
with respect to Hr and H^. This polynomial must identically vanish. Collecting similar terms 
we have 

iB{r)H + C{r) sin cl))h3„ - + (2/12 - B{r))h3 + C(r)(/i2 cos0 - sin<^) = 0, 
/13h« = 0, 2/l3,^ + h'^ - 2h'2h2 = 0. 

Thus, 

/i3(</,,if) = Q/i2 _ + /i4^ i/, /i'2 -/i^ -2cot<?!)/i2 + 2(l -/14) = 0, (33) 
where h^ is an arbitrary constant. 
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Clearly, that the Riccati equation (33) has the partial solution 

/i2 = tanc/) 

for /i4 = 1/2. 

Thus we find the solution of (31) 

h = H^^ + tan (j). 

The corresponding differential constraint h = has the general solution 

H = Hi{r)s\n(j) + H2{r), (34) 

where Hi and H2 are arbitrary functions. 

Next we use (34) and consider the equation (27) in more general form 

R^^ + r^Rrr + K{r)Rr + L{r)R + M(r) sin^ </> + N{r) sin (p + P(r) = 0, (35) 

where K{r), L{r), M{r), N{r), and P{r) are arbitrary functions. The BDEs method applied to 
the equation (35) gives rise to the following results: 

N{r) = 0, bi{r,(f)) = K{r), b2ir,<l)) = L{r) - 8 sm'^ 2(j), 
h = R(f,(j) — 2R^ cot 20. 

Integrating differential constraint h = corresponding to the BDE solution (34), we find 

R = Ri{r)sm'^ (f) + R2{r), (36) 

where Ri{r) and R2{r) are arbitrary functions. Using (34) and (36) we obtain the following 
corollary in terms of variables x, y, and z. 

Corollary 1. The following expressions for the unknown functions 



^2.-2 E I vj^_L_ 1 , = ( V- - 1 , (37) 



2 




Z 



2 



(Pi) = ^ , = ^^^^ + ^^^b-^ (38) 



flZ/oui lis to reduce the model (l)-(4) to a system of ODEs. 
Indeed, substituting (37) and (38) into (l)-(4) we obtain 
E'^E" E'^E' fa E' l\ 1 



2—--]-^{2{a-l)E + G-aTE')=Q, (39) 



G G \G E T J 

E^G" E^GWG' 4 1 ((2«-3)G + C.2$-arG')=0, (40) 



G G \G E T J G, y ' " S 



G G \G E T GpE"^) tG\G E 

E^R'l E'^R[ (G' ^E^_^_ _ E'^Ri /^2— - A— + '^^ ^'^ 



G G \G E T CipE'^J tG \ G E Gi^E^ 



Gp E'^H' / (H-l) 



Cip G \ T 



E^R'^ E'^R'^ fG' ^E' 1 a tG\ R2 f ^ G , ^ \ , ^E^Ri 



G G \G + G 

+ 2^(^-1)2^ = 0, (43) 

where r = y^^^ + rf. 
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R 0,6- 



(b) 
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Figure 1. Calculated profiles as ^ = 0: (a) normed profile of E, (b) normed profile of G, (c) profile 
of H, (d) normed profile of R. 



5 Calculation results 

The system of ODEs (39)-(43) has to satisfy the conditions 
E' = G' = H[ = R[ =R'^ = 0, T = 0, 



E — G — Hi — Ri — i?2 — 0, 



T 

T — )• OO. 



(44) 
(45) 



Conditions (44) take into account flow symmetry with respect to the OX axis. The boundary 
conditions (45) imply that all functions take zero values outside the turbulent wake. 

The system of ODEs (39)-(43) satisfying boundary conditions (44), (45) was solved numeri- 
cally. Additional difficulties are caused by the fact that the coefficients of ODEs have singula- 
rities. The problem was solved by the modified shooting method and asymptotic expansion of 
the solution in the vicinity of the singular point [3, 15] 



E = ci{t — a) r 



1 10/7 



G = -^^(T-a)i3/7 



7a 



+ o T - a 



13/7 



H 



7C, 



a{7C, - lOCe) 



(r - a) + o(|t - a|), 
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Figure 2. Calculated functions: (a) the function E/Eq, (b) the function G/Gq, (c) the function H, (d) 
the function R/Rq. 



Ri 



R2 



49C2(7(2a + l)Cp - 20aCe 



2a^{7Cp-WCey{5Ce-7Cp, 



+ T 



7C, 



2{5Ce 



i)2 + o(|r-a|2). 



The value of a is taken to be 0.23 in accordance with experimental data [9, 17]. The results 
for the problem solution are illustrated in Figs. 1 and 2. Fig. 1 shows the profiles of the functions 
E/Eq, G/Gq, H, and R/Rq as ^ = 0, where subscript denotes the axial value. The functions 
E/Eq, G/Gq, H, and R/Rq are plotted in Fig. 2. The functions E/Eq and G/Gq are bell- 
shaped and determine shapes of the normalized turbulent kinetic energy and the normalized 
kinetic energy dissipation rate respectively. Similarly, H and R/Rq determine shapes of the 
average density defect and the normalized density fluctuation varience respectively. 

The function H{0,ri) characterizing the degree of fluid mixing in the turbulent wake is given 

in Fig. Ic. The maximum value of this function equals 0.258. This is in consistent with the 

present notions of incomplete fluid mixing in the wakes [23, 24]. 

In Fig. 3 adapted from [6] the normalized values of the turbulent energy along the wake 
1/2 

axis Bq /Uq = e{x, 0, 0)/C/o are compared with experimental data [17], computational results [9] 
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Figure 3. Axial values of the turbulent energy. 

and results of numerical simulation based on two semi-empirical turbulence models (Modell and 
Model2 in [4, 6]). The coordinate x is normalized by the body diameter D. The results obtained 
here are in close agreement with Lin and Pao's experimental data. 
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